Data description

Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.

NB - if chain/sample parameters is modified - should be changed in the whole section(dataset)

Environment filtering 5 species

JSDM

data<-sim_data$EnvEvenSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

Y_cor<-cor(data$Y)

to_prec<-function(m){
  n<-dim(m)[1]
  Tau_n<-matrix(nrow=n, ncol=n)
  for (j in 1:n) {
   for (k in 1:n){
    Tau_n[j, k] <-  -m[j, k]/sqrt((m[j,j]*m[k,k]))
   }
  }
  return(Tau_n)
}

#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
cols = colorRampPalette(c("dark blue","white","red"))
col2 <- colorRampPalette(c("#4393C3", "#2166AC", "#053061",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#67001F", "#B2182B", "#D6604D", "#F4A582"))


corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")

#j_metric_e5<-metrics_jsdm(me5)
#cat("Success rate ")
#j_metric_e5$success_env
 cat(sprintf("Success rate: %s\n", metrics_jsdm(me5)$success_env))
## Success rate: 0.9

GJAM

##to setup chains parameters
data<-sim_data$EnvEvenSp5
#Rho_sign<-fit_gjam(data,2000,1000,"./gjam_models/gjam5env.rda")
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5env.rda")

#g_metric_e5<-metrics_gjam(Rho_sign)
#cat("Success rate ")
#g_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_gjam(Rho_sign)$success_env))
## Success rate: 0.8

HMSC

#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="Fit",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
  if (label=="Fit"){
    Y_data = subset(data, select = -env)
    ns<- ncol(Y_data)
    np <- nrow(Y_data)
    X<-scale(poly(data$env[1:np], 2))
    colnames(X)<-c("env","env2")
    studyDesign = data.frame(sample = as.factor(1:np))
    rL = HmscRandomLevel(units = studyDesign$sample)
    m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
             studyDesign = studyDesign, ranLevels = list(sample = rL))
    m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
    save(m, file = name)
    return(m)
  }
  if (label=="Load"){
    return(load_object(name))
  }
}




data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")

Convergence:

hm_conv<-function(mod){
  codaList = convertToCodaObject(mod)

  #convergence histograms
  hist(effectiveSize(codaList$Beta), main="ess(beta)",lwd=2,col=gray(.6))
  hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(beta)")
  
  hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)",lwd=2,col=gray(.6))
  hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf,lwd=2,col=gray(.6),      main="psrf(omega)")
}

hm_conv(hm_mod)

Study of interactions

metrics_hmsc<-function(Rho, comp=NULL,fac=NULL,only_env=T){
  #{
    if(!is.null(comp) & !is.null(fac)){
    fac_comp <- fac-comp
    TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
    for(i in 1:(nrow(Rho)-1)){
      for(j in (i+1):ncol(Rho)){
        if(fac_comp[i,j]>0 & Rho[i,j] >0)
          TP_fac=TP_fac+1
        if(fac_comp[i,j]>0 & Rho[i,j] == 0)
          FN_fac=FN_fac+1
        if(fac_comp[i,j]==0 & Rho[i,j] == 0)
           TN=TN+1
        if(fac_comp[i,j]==0 & Rho[i,j] != 0)
           FP=FP+1
        if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
          TP_comp=TP_comp+1
        if(fac_comp[i,j]<0 & Rho[i,j] == 0)
          FN_comp=FN_comp+1
        if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
          wrong=wrong+1
      }
    }
    
  }else{
    if(!is.null(comp)){
    TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
    TP_fac <- FN_fac<-NULL
    for(i in 1:(nrow(Rho)-1)){
      for(j in (i+1):ncol(Rho)){
        if(comp[i,j]>0 & Rho[i,j] <0)
          TP_comp=TP_comp+1
        if(comp[i,j]>0 & Rho[i,j] >0)
          wrong=wrong+1
        if(comp[i,j]>0 & Rho[i,j] == 0)
          FN_comp=FN_comp+1
        if(comp[i,j]==0 & Rho[i,j] == 0)
          TN=TN+1
        if(comp[i,j]==0 & Rho[i,j] != 0)
          FP=FP+1
      }
    }
  #  }
  #  }
  }
  if(!is.null(fac)){
    
            TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
            TP_comp <- FN_comp<-NULL
            for(i in 1:(nrow(Rho)-1)){
              for(j in (i+1):ncol(Rho)){
                   if(fac[i,j]>0 & Rho[i,j] >0)
                      TP_fac=TP_fac+1
                   if(fac[i,j]>0 & Rho[i,j] <0)
                      wrong = wrong+1
                   if(fac[i,j]>0 & Rho[i,j] == 0)
                      FN_fac=FN_fac+1
                   if(fac[i,j]==0 & Rho[i,j] == 0)
                      TN=TN+1
                   if(fac[i,j]==0 & Rho[i,j] != 0)
                      FP=FP+1
              }
            }
  }
  }  
  
  if(only_env){
    
      FP <- TN <- 0
      TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
         for(i in 1:(nrow(Rho)-1)){
            for(j in (i+1):nrow(Rho)){
                if(Rho[i,j] != 0)
                      FP=FP+1
                if(Rho[i,j] == 0)
                      TN=TN+1
              }
            }
      
  }
  success_env<-TN/(TN+FP)
  
  if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
  
  if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
  
  list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
       "wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}



hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
  getOmega = function(a,r=1)
  return(crossprod(a$Lambda[[r]]))
  ns<-mod$ns
  postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
  postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
  
  postOmega<-abind(postOmega1,postOmega2,along=3)
  postOmegaMean = apply(postOmega,c(1,2),mean)
  postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
  postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
  
  postR<-array(dim=c(ns,ns,nchains*nsamples))
  for(i in 1:dim(postOmega)[3])
  postR[,,i]<-stats::cov2cor(postOmega[,,i])
  
  
  postRMean = apply(postR,c(1,2),mean)
  postRUp=apply(postR,c(1,2),quantile,0.95)
  postRLo=apply(postR,c(1,2),quantile,0.05)
  
  Tau = solve(postOmegaMean)
  Tau_n = -cov2cor(Tau)
  
  
  Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
  
  # Omegacor<- computeAssociations(m)
  # supportLevel<- 0.95
  # toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
  # corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
  
  par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  corrplot(cor(hm_mod$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(postRMean, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R")
  corrplot(Toplot_R, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Plot only non zero value")
  corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation matrix")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
  
  Toplot_R
}

Rho_sign<-hm_inter(hm_mod)

# h_metric_e5<-metrics_hmsc(Rho_sign)
# cat("Success rate ")
# h_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_hmsc(Rho_sign)$success_env))
## Success rate: 1

Environmental filtering 10 species

JSDM

me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)
## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
jsdm_conv(me10)
## Maximum Rhat value for Beta: 1.5602
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.7953

me10$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
j_metric_e10<-metrics_jsdm(me10)
cat("Success rate ")
## Success rate
j_metric_e10$success_env
## [1] 0.9333333
data<-sim_data$EnvEvenSp10
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
  par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
  corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("EnvRho")
  corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("EnvRho signif")
  corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Rho")
  corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Rho signif")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}

plot_cor_jsdm(me10,data$Y)

Gjam

data<-sim_data$EnvEvenSp10

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10env.rda")
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10env.rda")

#gje10<-load_object("./gjam_models/gjam10env.rda")

g_metric_e10<-metrics_gjam(Rho_sign)
cat("Success rate ")
## Success rate
g_metric_e10$success_env
## [1] 0.2666667
#to check posterior density of s in Sigma 
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))

HMSC

data<-sim_data$EnvEvenSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10env.rda" )

hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod)

h_metric_e10<-metrics_hmsc(Rho_sign)
# cat("Success rate ")
# h_metric_e10$success_env
# h_metric_e10$success_fac

cat(" ")
cat(sprintf("Success rate: %s\n", h_metric_e10$success_env))
## Success rate: 1
cat(sprintf("Success rate: %s\n", h_metric_e10$success_fac))

Environmental filtering 20 species

JSDM

me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)
## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
me20$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
jsdm_conv(me20)
## Maximum Rhat value for Beta: 1.2174
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3236

j_metric_e20<-metrics_jsdm(me20)
cat("Success rate ")
## Success rate
j_metric_e20$success_env
## [1] 0.9894737
data<-sim_data$EnvEvenSp20
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################

plot_cor_jsdm(me20,data$Y)

Gjam

data<-sim_data$EnvEvenSp20

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

Rho_sign<-load_gjam(data,name="./gjam_models/gjam20env.rda")

#gje20<-load_object("./gjam_models/gjam20env.rda")

# g_metric_e20<-metrics_gjam(Rho_sign)
# cat("Success rate ")
# g_metric_e20$success_env

cat(" ")
cat(sprintf("Success rate env: %s\n", round(metrics_gjam(Rho_sign)$success_env,3)))
## Success rate env: 0.384
#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))

Gjam dimension reduction

data<-sim_data$EnvEvenSp20

#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

#load_gjam(data,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")


#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
 data <- list(
    Y = subset(data, select = -env),
    X = cbind(1, scale(poly(data$env, 2))),
    covx = cov(cbind(1, scale(poly(data$env, 2)))),
    K = 3,
    J = ncol(data) - 1,
    n = nrow(data),
    I = diag(ncol(data) - 1),
    df = ncol(data)
  )
xdata<-as.data.frame(data$X[,-1])
colnames(xdata)<- c("env","env2")
ydata<-as.data.frame(data$Y)
###formula
rl   <- list(r = 8, N = 20)
formula<-as.formula( ~env+ env2)
ml  <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
####fit


mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
## 
## Observations and responses:
## [1] 500  20
## Warning in .setupReduct(modelList, S, Q, n): dimension reduction requires
## reductList$N < no. responses
## 
## Dimension reduced from 20 X 20 -> 20 X 8 responses
## ===========================================================================
## expanding covariance chains
## ===========================================================================
save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
summary(mod_gjam1)
## 
## Sensitivity by predictor variables f:
##      Estimate     SE CI_025 CI_975
## env     0.958 0.0836  0.823   1.16
## env2    1.310 0.0685  1.180   1.45
## 
## Coefficient matrix B:
##             sp01   sp02    sp03   sp04   sp05   sp06   sp07   sp08   sp09
## intercept -0.341 -0.101  0.0155  0.114  0.292  0.337  0.470  0.482  0.484
## env       -0.473 -0.480 -0.4800 -0.474 -0.470 -0.460 -0.437 -0.383 -0.286
## env2       0.285  0.106  0.0122 -0.104 -0.255 -0.325 -0.438 -0.455 -0.478
## RMSPE      0.406  0.405  0.4090  0.411  0.404  0.403  0.409  0.423  0.434
##              sp10   sp11   sp12   sp13   sp14   sp15   sp16   sp17
## intercept  0.4880  0.488  0.486  0.478  0.476  0.395  0.272  0.175
## env       -0.0399  0.207  0.457  0.437  0.453  0.466  0.477  0.478
## env2      -0.4830 -0.474 -0.481 -0.478 -0.457 -0.354 -0.262 -0.175
## RMSPE      0.4430  0.446  0.428  0.425  0.405  0.407  0.407  0.409
##               sp18   sp19   sp20
## intercept -0.00278 -0.092 -0.276
## env        0.48000  0.479  0.470
## env2       0.00499  0.092  0.237
## RMSPE      0.40400  0.407  0.404
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept -0.34100 0.0696 -0.47300 -0.2060     *
## sp01_env       -0.47300 0.0268 -0.49900 -0.3980     *
## sp01_env2       0.28500 0.0683  0.14900  0.4210     *
## sp02_intercept -0.10100 0.0722 -0.23800  0.0340      
## sp02_env       -0.48000 0.0235 -0.50000 -0.4170     *
## sp02_env2       0.10600 0.0674 -0.02310  0.2410      
## sp03_intercept  0.01550 0.0734 -0.12500  0.1600      
## sp03_env       -0.48000 0.0210 -0.49900 -0.4240     *
## sp03_env2       0.01220 0.0677 -0.11900  0.1410      
## sp04_intercept  0.11400 0.0751 -0.02790  0.2590      
## sp04_env       -0.47400 0.0306 -0.49900 -0.3990     *
## sp04_env2      -0.10400 0.0677 -0.23200  0.0286      
## sp05_intercept  0.29200 0.0767  0.14100  0.4440     *
## sp05_env       -0.47000 0.0288 -0.49900 -0.3930     *
## sp05_env2      -0.25500 0.0710 -0.39300 -0.1130     *
## sp06_intercept  0.33700 0.0692  0.19200  0.4660     *
## sp06_env       -0.46000 0.0387 -0.49900 -0.3620     *
## sp06_env2      -0.32500 0.0692 -0.46300 -0.1900     *
## sp07_intercept  0.47000 0.0270  0.40000  0.4990     *
## sp07_env       -0.43700 0.0667 -0.49900 -0.2660     *
## sp07_env2      -0.43800 0.0483 -0.49700 -0.3220     *
## sp08_intercept  0.48200 0.0176  0.43300  0.5000     *
## sp08_env       -0.38300 0.0949 -0.49500 -0.1350     *
## sp08_env2      -0.45500 0.0369 -0.49900 -0.3630     *
## sp09_intercept  0.48400 0.0155  0.44200  0.5000     *
## sp09_env       -0.28600 0.1140 -0.48700 -0.0525     *
## sp09_env2      -0.47800 0.0208 -0.49900 -0.4220     *
## sp10_intercept  0.48800 0.0116  0.45800  0.5000     *
## sp10_env       -0.03990 0.1210 -0.25200  0.2070      
## sp10_env2      -0.48300 0.0165 -0.50000 -0.4370     *
## sp11_intercept  0.48800 0.0116  0.45600  0.5000     *
## sp11_env        0.20700 0.1130  0.00925  0.4570     *
## sp11_env2      -0.47400 0.0241 -0.49900 -0.4070     *
## sp12_intercept  0.48600 0.0139  0.44800  0.5000     *
## sp12_env        0.45700 0.0400  0.34600  0.4990     *
## sp12_env2      -0.48100 0.0178 -0.50000 -0.4330     *
## sp13_intercept  0.47800 0.0210  0.42100  0.4990     *
## sp13_env        0.43700 0.0596  0.28400  0.4980     *
## sp13_env2      -0.47800 0.0214 -0.49900 -0.4210     *
## sp14_intercept  0.47600 0.0238  0.41300  0.4990     *
## sp14_env        0.45300 0.0529  0.31300  0.4990     *
## sp14_env2      -0.45700 0.0356 -0.49800 -0.3710     *
## sp15_intercept  0.39500 0.0608  0.26300  0.4920     *
## sp15_env        0.46600 0.0318  0.38100  0.4990     *
## sp15_env2      -0.35400 0.0682 -0.47700 -0.2120     *
## sp16_intercept  0.27200 0.0747  0.12700  0.4160     *
## sp16_env        0.47700 0.0259  0.40900  0.4990     *
## sp16_env2      -0.26200 0.0734 -0.40900 -0.1180     *
## sp17_intercept  0.17500 0.0756  0.02290  0.3210     *
## sp17_env        0.47800 0.0228  0.41300  0.4990     *
## sp17_env2      -0.17500 0.0720 -0.31000 -0.0364     *
## sp18_intercept -0.00278 0.0748 -0.14300  0.1420      
## sp18_env        0.48000 0.0226  0.42700  0.4990     *
## sp18_env2       0.00499 0.0703 -0.13200  0.1430      
## sp19_intercept -0.09200 0.0729 -0.23500  0.0550      
## sp19_env        0.47900 0.0231  0.42000  0.4990     *
## sp19_env2       0.09200 0.0693 -0.04010  0.2280      
## sp20_intercept -0.27600 0.0756 -0.42900 -0.1320     *
## sp20_env        0.47000 0.0320  0.39100  0.4990     *
## sp20_env2       0.23700 0.0698  0.10500  0.3780     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.415, and the DIC is 77360.  Computation involved 2500 Gibbs steps, with a burnin of 500.  Dimension reduction was implemented with N = 11 and r = 8.
Tau <- solve(mod_gjam1$parameters$sigMu)
Tau_n = to_prec(Tau)

#postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
#postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)


#pH<-convert_to_m(postH)
#pL<-convert_to_m(postL)

#R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))

g_metric_e20_dim_red<-metrics_gjam(Rho_sign)
cat("Success rate ")
## Success rate
g_metric_e20_dim_red$success_env
## [1] 0.3842105
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("E matrix")
# corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("Tau")
# corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("R signif")
corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")

data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod)

# h_metric_e20<-metrics_hmsc(Rho_sign)
# cat("Success score ")
# h_metric_e20$success_env
cat(" ")
cat(sprintf("Success rate non-interacting: %s\n", round(metrics_hmsc(Rho_sign)$success_env,3)))
## Success rate non-interacting: 1

Environmental filtering + Facilitation dense 5 species

JSDM

mfd5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mfd5)
## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#mfd5$Rhat
jsdm_conv(mfd5)
## Maximum Rhat value for Beta: 1.3081
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1421

mfd5$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
j_metric_facDense5<-metrics_jsdm(mfd5,fac = fac_inter[[4]],only_env = F)
cat("Success rate ")
## Success rate
j_metric_facDense5$success_env
## [1] 1
j_metric_facDense5$success_fac
## [1] 0
data<-sim_data$FacDenseSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm(mfd5,data$Y,fac_inter[[4]])

Gjam

data<-sim_data$FacDenseSp5
#Rho_sign<-fit_gjam(data,10000,2000,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])

#gjfd5<-load_object("./gjam_models/gjam5f.rda")

g_metric_facDense5<-metrics_gjam(Rho_sign,fac=fac_inter[[4]], only_env = F)
#cat("Success rate ")
#g_metric_facDense5$success_env
#g_metric_facDense5$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 0.375
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0.5

HMSC

data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[4]])
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1
Rho_sign<-hm_inter(hm_mod)

h_metric_facDense5<-metrics_hmsc(Rho_sign,fac=fac_inter[[4]],only_env = F)
#cat("Success rate ")
# h_metric_facDense5$success_env
# h_metric_facDense5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Facilitation dense 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 0.925
## [1] 0

data<-sim_data$FacDenseSp10

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
load_gjam(data,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])

##            sp01       sp02       sp03       sp04       sp05       sp06
## sp01  1.0000000 -0.3904084  0.0000000  0.2493194 -0.1879177 -0.6181350
## sp02 -0.3904084  1.0000000  0.0000000 -0.3751301  0.0000000  0.0000000
## sp03  0.0000000  0.0000000  1.0000000  0.2974028  0.0000000  0.0000000
## sp04  0.2493194 -0.3751301  0.2974028  1.0000000  0.0000000  0.3768348
## sp05 -0.1879177  0.0000000  0.0000000  0.0000000  1.0000000  0.0000000
## sp06 -0.6181350  0.0000000  0.0000000  0.3768348  0.0000000  1.0000000
## sp07  0.0000000  0.0000000  0.0000000 -0.7225666 -0.1509510 -0.3458261
## sp08  0.0000000  0.4069187  0.0000000  0.0000000  0.0000000  0.0000000
## sp09 -0.1450869  0.0000000 -0.1989728  0.0000000  0.0000000  0.0000000
## sp10 -0.5235049 -0.2648058  0.0000000  0.0000000  0.0000000  0.6694006
##            sp07      sp08       sp09       sp10
## sp01  0.0000000 0.0000000 -0.1450869 -0.5235049
## sp02  0.0000000 0.4069187  0.0000000 -0.2648058
## sp03  0.0000000 0.0000000 -0.1989728  0.0000000
## sp04 -0.7225666 0.0000000  0.0000000  0.0000000
## sp05 -0.1509510 0.0000000  0.0000000  0.0000000
## sp06 -0.3458261 0.0000000  0.0000000  0.6694006
## sp07  1.0000000 0.0000000  0.3544135  0.3680399
## sp08  0.0000000 1.0000000  0.0000000  0.0000000
## sp09  0.3544135 0.0000000  1.0000000  0.3121062
## sp10  0.3680399 0.0000000  0.3121062  1.0000000
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")

g_metric_facDense10<-metrics_gjam(Rho_sign,fac=fac_inter[[5]], only_env = F)
# cat("Success rate ")
# g_metric_facDense10$success_env
# g_metric_facDense10$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense10$success_fac,3)))
## Success rate for facilitation: 0

HMSC

## Success rate for non-interacting: 1
## Success rate for facilitation: 0

Environmental filtering + Facilitation dense 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9833333
## [1] 0

Gjam

data<-sim_data$FacDenseSp20

#Rho_sign<-fit_gjam(data,10000,1500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])

#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

g_metric_facDense20<-metrics_gjam(Rho_sign,fac=fac_inter[[6]], only_env = F)
# cat("Success rate ")
# g_metric_facDense20$success_env
# g_metric_facDense20$success_fac
# 
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 0.411
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0.5

HMSC

data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, interact = fac_inter[[6]])

h_metric_facDense20<-metrics_hmsc(Rho_sign,fac=fac_inter[[6]],only_env = F)
# cat("Success rate ")
# h_metric_facDense20$success_env
# h_metric_facDense20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Facilitation sparse 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.889
## Success rate for facilitation: 1

Gjam

data<-sim_data$FacSparseSp5

#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam5fs.rda",interact=fac_inter[[7]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5fs.rda", interact=fac_inter[[7]])

#gjfd5<-load_object("./gjam_models/gjam5fs.rda")

g_metric_facSparse5<-metrics_gjam(Rho_sign,fac=fac_inter[[7]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse5$success_env
# g_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacSparseSp5

hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fs.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, interact = fac_inter[[7]])
h_metric_facSparse5<-metrics_hmsc(Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Facilitation sparse 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.93
## Success rate for facilitation: 0

data<-sim_data$FacSparseSp10

#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam10fs.rda",interact=fac_inter[[8]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10fs.rda", interact=fac_inter[[8]])

#gjfd5<-load_object("./gjam_models/gjam10fs.rda")

g_metric_facSparse10<-metrics_gjam(Rho_sign,fac=fac_inter[[8]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse10$success_env
# g_metric_facSparse10$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 0.372
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0

HMSC

data<-sim_data$FacSparseSp10

hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10fs.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, interact = fac_inter[[8]])
h_metric_facSparse10<-metrics_hmsc(Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse10$success_env
# h_metric_facSparse10$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Facilitation sparse 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.978
## Success rate for facilitation: 0

Gjam

data<-sim_data$FacSparseSp20

#fit_gjam(data,10000,1500,"./gjam_models/gjam20fs.rda",interact=fac_inter[[9]])
load_gjam(data,name="./gjam_models/gjam20fs.rda", interact=fac_inter[[9]])

##            sp01       sp02       sp03       sp04       sp05       sp06
## sp01  1.0000000  0.0000000  0.0000000  0.2481647 -0.1405985 -0.1778950
## sp02  0.0000000  1.0000000 -0.1750012  0.0000000  0.0000000 -0.6945088
## sp03  0.0000000 -0.1750012  1.0000000  0.0000000  0.1633683  0.3398000
## sp04  0.2481647  0.0000000  0.0000000  1.0000000  0.1812365  0.0000000
## sp05 -0.1405985  0.0000000  0.1633683  0.1812365  1.0000000  0.0000000
## sp06 -0.1778950 -0.6945088  0.3398000  0.0000000  0.0000000  1.0000000
## sp07 -0.2262923  0.0000000  0.2351144 -0.2474855 -0.3272462  0.0000000
## sp08  0.0000000  0.4861375  0.0000000 -0.2295886  0.0000000 -0.3307800
## sp09  0.0000000  0.0000000  0.5327952  0.0000000  0.0000000  0.3040397
## sp10  0.2989764  0.0000000  0.4502637  0.0000000  0.0000000 -0.1168824
## sp11  0.3747704 -0.2002289  0.3865259  0.2438452 -0.2970431  0.0000000
## sp12 -0.2809890  0.0000000  0.1516279 -0.4507748  0.0000000 -0.1320808
## sp13  0.0000000 -0.2293087  0.3528826  0.0000000  0.2622547  0.2939977
## sp14  0.0000000  0.1499564  0.0000000  0.0000000  0.0000000  0.0000000
## sp15  0.3873898  0.0000000 -0.2520398 -0.2776349  0.1619795  0.0000000
## sp16 -0.2582433  0.0000000 -0.1597638 -0.2483336 -0.2604524  0.2649414
## sp17  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.2525415
## sp18  0.0000000 -0.1554376  0.0000000  0.0000000  0.0000000  0.0000000
## sp19 -0.1567031  0.2745205  0.1115394  0.0000000  0.0000000  0.0000000
## sp20 -0.3296659  0.0000000  0.0000000  0.0000000  0.1530916 -0.2474593
##            sp07       sp08       sp09       sp10       sp11       sp12
## sp01 -0.2262923  0.0000000  0.0000000  0.2989764  0.3747704 -0.2809890
## sp02  0.0000000  0.4861375  0.0000000  0.0000000 -0.2002289  0.0000000
## sp03  0.2351144  0.0000000  0.5327952  0.4502637  0.3865259  0.1516279
## sp04 -0.2474855 -0.2295886  0.0000000  0.0000000  0.2438452 -0.4507748
## sp05 -0.3272462  0.0000000  0.0000000  0.0000000 -0.2970431  0.0000000
## sp06  0.0000000 -0.3307800  0.3040397 -0.1168824  0.0000000 -0.1320808
## sp07  1.0000000  0.1866792  0.1815724 -0.1164417  0.0000000 -0.1966267
## sp08  0.1866792  1.0000000  0.2297215  0.0000000  0.0000000  0.0000000
## sp09  0.1815724  0.2297215  1.0000000  0.1862230  0.0000000  0.0000000
## sp10 -0.1164417  0.0000000  0.1862230  1.0000000  0.1139547  0.0000000
## sp11  0.0000000  0.0000000  0.0000000  0.1139547  1.0000000  0.0000000
## sp12 -0.1966267  0.0000000  0.0000000  0.0000000  0.0000000  1.0000000
## sp13  0.0000000 -0.4251175 -0.1844692 -0.1109607  0.0000000  0.0000000
## sp14  0.0000000 -0.3169754  0.0000000 -0.1806383  0.1229519  0.0000000
## sp15 -0.2393655  0.3940648  0.0000000  0.0000000 -0.1572325  0.0000000
## sp16  0.0000000  0.0000000 -0.2478583 -0.1822943  0.0000000  0.0000000
## sp17 -0.4179798 -0.2271328  0.2818370 -0.3565917 -0.2899321  0.0000000
## sp18 -0.3452972  0.0000000  0.3590445 -0.1968948  0.0000000  0.0000000
## sp19  0.0000000 -0.3549123 -0.3927233  0.0000000  0.0000000  0.2486426
## sp20 -0.2687936  0.0000000  0.0000000 -0.1423872  0.3274940  0.4736903
##            sp13       sp14       sp15       sp16       sp17       sp18
## sp01  0.0000000  0.0000000  0.3873898 -0.2582433  0.0000000  0.0000000
## sp02 -0.2293087  0.1499564  0.0000000  0.0000000  0.0000000 -0.1554376
## sp03  0.3528826  0.0000000 -0.2520398 -0.1597638  0.0000000  0.0000000
## sp04  0.0000000  0.0000000 -0.2776349 -0.2483336  0.0000000  0.0000000
## sp05  0.2622547  0.0000000  0.1619795 -0.2604524  0.0000000  0.0000000
## sp06  0.2939977  0.0000000  0.0000000  0.2649414  0.2525415  0.0000000
## sp07  0.0000000  0.0000000 -0.2393655  0.0000000 -0.4179798 -0.3452972
## sp08 -0.4251175 -0.3169754  0.3940648  0.0000000 -0.2271328  0.0000000
## sp09 -0.1844692  0.0000000  0.0000000 -0.2478583  0.2818370  0.3590445
## sp10 -0.1109607 -0.1806383  0.0000000 -0.1822943 -0.3565917 -0.1968948
## sp11  0.0000000  0.1229519 -0.1572325  0.0000000 -0.2899321  0.0000000
## sp12  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
## sp13  1.0000000  0.1641325  0.0000000  0.0000000  0.0000000 -0.1065072
## sp14  0.1641325  1.0000000  0.0000000  0.2800982  0.2165972  0.0000000
## sp15  0.0000000  0.0000000  1.0000000  0.0000000  0.2232036  0.0000000
## sp16  0.0000000  0.2800982  0.0000000  1.0000000  0.0000000  0.2704149
## sp17  0.0000000  0.2165972  0.2232036  0.0000000  1.0000000  0.2673245
## sp18 -0.1065072  0.0000000  0.0000000  0.2704149  0.2673245  1.0000000
## sp19  0.6487585  0.0000000 -0.2207364  0.0000000  0.0000000 -0.3499635
## sp20 -0.2588774 -0.4217491 -0.2976427  0.0000000  0.0000000  0.2694292
##            sp19       sp20
## sp01 -0.1567031 -0.3296659
## sp02  0.2745205  0.0000000
## sp03  0.1115394  0.0000000
## sp04  0.0000000  0.0000000
## sp05  0.0000000  0.1530916
## sp06  0.0000000 -0.2474593
## sp07  0.0000000 -0.2687936
## sp08 -0.3549123  0.0000000
## sp09 -0.3927233  0.0000000
## sp10  0.0000000 -0.1423872
## sp11  0.0000000  0.3274940
## sp12  0.2486426  0.4736903
## sp13  0.6487585 -0.2588774
## sp14  0.0000000 -0.4217491
## sp15 -0.2207364 -0.2976427
## sp16  0.0000000  0.0000000
## sp17  0.0000000  0.0000000
## sp18 -0.3499635  0.2694292
## sp19  1.0000000  0.0000000
## sp20  0.0000000  1.0000000
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

g_metric_facSparse20<-metrics_gjam(Rho_sign,fac=fac_inter[[9]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse20$success_env
# g_metric_facSparse20$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0

HMSC

HMSC

data<-sim_data$FacSparseSp20
hm_mod<-fit_hmsc(data,nsamples=3000, nchains=2,"Load",name="./HMmodels/hm20fs.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=3000, nchains=2,interact = fac_inter[[9]])
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1    0    0    0    0    0    0    0    0     0     0     0     0
##  [2,]    0    1    0    0    0    0    0    0    0     0     0     0     0
##  [3,]    0    0    1    0    0    0    0    0    0     0     0     0     0
##  [4,]    0    0    0    1    0    0    0    0    0     0     0     0     0
##  [5,]    0    0    0    0    1    0    0    0    0     0     0     0     0
##  [6,]    0    0    0    0    0    1    0    0    0     0     0     0     0
##  [7,]    0    0    0    0    0    0    1    0    0     0     0     0     0
##  [8,]    0    0    0    0    0    0    0    1    0     0     0     0     0
##  [9,]    0    0    0    0    0    0    0    0    1     0     0     0     0
## [10,]    0    0    0    0    0    0    0    0    0     1     0     0     0
## [11,]    0    0    0    0    0    0    0    0    0     0     1     0     0
## [12,]    0    0    0    0    0    0    0    0    0     0     0     1     0
## [13,]    0    0    0    0    0    0    0    0    0     0     0     0     1
## [14,]    0    0    0    0    0    0    0    0    0     0     0     0     0
## [15,]    0    0    0    0    0    0    0    0    0     0     0     0     0
## [16,]    0    0    0    0    0    0    0    0    0     0     0     0     0
## [17,]    0    0    0    0    0    0    0    0    0     0     0     0     0
## [18,]    0    0    0    0    0    0    0    0    0     0     0     0     0
## [19,]    0    0    0    0    0    0    0    0    0     0     0     0     0
## [20,]    0    0    0    0    0    0    0    0    0     0     0     0     0
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]     0     0     0     0     0     0     0
##  [2,]     0     0     0     0     0     0     0
##  [3,]     0     0     0     0     0     0     0
##  [4,]     0     0     0     0     0     0     0
##  [5,]     0     0     0     0     0     0     0
##  [6,]     0     0     0     0     0     0     0
##  [7,]     0     0     0     0     0     0     0
##  [8,]     0     0     0     0     0     0     0
##  [9,]     0     0     0     0     0     0     0
## [10,]     0     0     0     0     0     0     0
## [11,]     0     0     0     0     0     0     0
## [12,]     0     0     0     0     0     0     0
## [13,]     0     0     0     0     0     0     0
## [14,]     1     0     0     0     0     0     0
## [15,]     0     1     0     0     0     0     0
## [16,]     0     0     1     0     0     0     0
## [17,]     0     0     0     1     0     0     0
## [18,]     0     0     0     0     1     0     0
## [19,]     0     0     0     0     0     1     0
## [20,]     0     0     0     0     0     0     1
h_metric_facSparse20<-metrics_hmsc(Rho_sign,fac=fac_inter[[9]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse20$success_env
# h_metric_facSparse20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Competition dense 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## Successful convergence based on Rhat values (all < 1.1).
## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 1
## [1] 1
## Success rate for non-interacting: 1
## Success rate for facilitation: 0

Gjam

data<-sim_data$CompDenseSp5

#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam5cmpd.rda",interact=comp_inter[[10]])

Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmpd.rda", interact=(-1)*comp_inter[[10]])

#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")

g_metric_compDense5<-metrics_gjam(Rho_sign,comp=comp_inter[[10]], only_env = F)
# cat("Success rate ")
# g_metric_compDense5$success_env
# g_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 0.625
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompDenseSp5

hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hm5cmpd.rda" )
#hm_mod<-fit_hmsc(data,"Fit",nsamples=3000, nchains=2 )

hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[10]])

h_metric_compDense5<-metrics_hmsc(Rho_sign,comp=comp_inter[[10]],only_env = F)
# cat("Success rate ")
# h_metric_compDense5$success_env
# h_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1

Environmental filtering + Competition dense 10 species

## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.875
## Success rate for facilitation: 1

Gjam

data<-sim_data$CompDenseSp10

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10cmpd.rda",interact=comp_inter[[11]])

Rho_sign<-load_gjam(data,name="./gjam_models/gjam10cmpd.rda", interact= (-1)*comp_inter[[11]])

#gjfd5<-load_object("./gjam_models/gjam10cmpd.rda")

g_metric_compDense10<-metrics_gjam(Rho_sign,comp=comp_inter[[11]], only_env = F)
# cat("Success rate ")
# g_metric_compDense10$success_env
# g_metric_compDense10$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.35
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompDenseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm10cmpd.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod,nsamples=5000, nchains=2, interact =  (-1)*comp_inter[[11]])

h_metric_compDense10<-metrics_hmsc(Rho_sign,comp=comp_inter[[11]],only_env = F)
# cat("Success rate ")
# h_metric_compDense10$success_env
# h_metric_compDense10$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.875
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense10$success_comp,3)))
## Success rate for competition: 1

Environmental filtering + Competition dense 20 species

## Summary for model '/tmp/RtmpKix4lZ/file40e96843124a'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2061.424 minutes at time 2019-04-16 03:39:17.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.9065
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.6944

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.956
## Success rate for competition: 0.7

Gjam

data<-sim_data$CompDenseSp20

#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam20cmpd.rda",interact= (-1)*comp_inter[[12]])

Rho_sign<-load_gjam(data,name="./gjam_models/gjam20cmpd.rda", interact= (-1)*comp_inter[[12]])

#gjfd5<-load_object("./gjam_models/gjam20cmpd.rda")

g_metric_compDense20<-metrics_gjam(Rho_sign,comp=comp_inter[[12]], only_env = F)
# cat("Success rate ")
# g_metric_compDense20$success_env
# g_metric_compDense20$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.406
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.9
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompDenseSp20

hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm20cmpd.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=5000, nchains=2,interact =  (-1)*comp_inter[[12]])

h_metric_compDense20<-metrics_hmsc(Rho_sign,comp=comp_inter[[12]],only_env = F)
# cat("Success rate ")
# h_metric_compDense20$success_env
# h_metric_compDense20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.967
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.4

Environmental filtering + Competition sparse 5 species

## Summary for model '/tmp/RtmpKix4lZ/file40e97e225117'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.251 minutes at time 2019-04-17 14:00:46.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1876
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2038

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 1
## Success rate for competition: 1

Gjam

data<-sim_data$CompSparseSp5


#Rho_sign<-fit_gjam(data,5000,2000,"./gjam_models/gjam5cmps.rda",interact= (-1)*comp_inter[[13]])

Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmps.rda", interact= (-1)*comp_inter[[13]])

#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")

g_metric_compSparse5<-metrics_gjam(Rho_sign,comp=comp_inter[[13]], only_env = F)
cat("Success rate ")
## Success rate
g_metric_compSparse5$success_env
## [1] 0.7777778
g_metric_compSparse5$success_comp
## [1] 1
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse5$success_env,3)))
## Success rate for non-interacting: 0.778
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse5$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompSparseSp5

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm5cmps.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[13]])

h_metric_facSparse5<-metrics_hmsc(Rho_sign,comp=comp_inter[[13]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_facSparse5$success_comp,3)))
## Success rate for competition: 1

Environmental filtering + Competition sparse 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e978d641ed'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 685.719 minutes at time 2019-04-17 14:30:01.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.116
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1753

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9534884
## [1] 1
## Success rate for non-interacting: 0.953
## Success rate for competition: 1

Gjam

data<-sim_data$CompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam10cmps.rda",interact= (-1)*comp_inter[[14]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10cmps.rda", interact= (-1)*comp_inter[[14]])

g_metric_compSparse10<-metrics_gjam(Rho_sign,comp=comp_inter[[14]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse10$success_env
# g_metric_compSparse10$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse10$success_env,3)))
## Success rate for non-interacting: 0.349
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm10cmps.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[14]])

h_metric_facSparse10<-metrics_hmsc(Rho_sign,comp=comp_inter[[14]],only_env = F)
cat("Success rate ")
## Success rate
h_metric_facSparse10$success_env
## [1] 0.9534884
h_metric_facSparse10$success_comp
## [1] 0.5
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 0.953
cat(sprintf("Success rate for competition: %s\n", round(h_metric_facSparse10$success_comp,3)))
## Success rate for competition: 0.5

Environmental filtering + Competition sparse 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9379ee963'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2084.707 minutes at time 2019-04-18 01:55:46.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2652
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3883

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.968
## Success rate for competition: 0

Gjam

data<-sim_data$CompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam20cmps.rda",interact= (-1)*comp_inter[[15]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam20cmps.rda", interact= (-1)*comp_inter[[15]])

g_metric_compSparse20<-metrics_gjam(Rho_sign,comp=comp_inter[[15]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse20$success_env
# g_metric_compSparse20$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 0.296
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0.75
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompSparseSp20

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm20cmps.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[15]])

h_metric_compSparse20<-metrics_hmsc(Rho_sign,comp=comp_inter[[15]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0

Environmental filtering + Competition +Facilitation 5 dense species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e93b49ad29'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 28.267 minutes at time 2019-04-19 12:40:31.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1667
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2086

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.5
## Success rate for facilitation: 0.5

Gjam

data<-sim_data$FacCompDenseSp5
#Rho_sign<-fit_gjam(data,30000,15000,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
Rho_sign2<-load_gjam(data,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])

#Rho_sign3<-fit_gjam(data,5000,1000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign4<-fit_gjam(data,5000,2500,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign5<-fit_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmp_facd.rda", (-1)*comp_inter[[16]]+fac_inter[[16]])

g_metric_cmp_facd5<-metrics_gjam(Rho_sign2,comp=comp_inter[[16]], fac=fac_inter[[16]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_cmp_facd5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(g_metric_cmp_facd5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_cmp_facd5$success_fac,3)))
## Success rate for facilitation: 1
# par(mfrow=c(2,3))
#corrplot((Rho_sign+Rho_sign2)/2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign3, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign4, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign5, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")


#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacCompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hmcmp_facd5.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact =  (-1)*comp_inter[[16]]+ fac_inter[[16]])

h_metric_FacCompDenseSp5<-metrics_hmsc(Rho_sign,comp=comp_inter[[16]],fac= fac_inter[[16]],only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDenseSp5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDenseSp5$success_comp,3)))
## Success rate for competition: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDenseSp5$success_fac,3)))
## Success rate for facilitation: 0.5

Environmental filtering + Competition +Facilitation 10 dense species

JSDM

Gjam

data<-sim_data$FacCompDenseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd10.rda",interact=  (-1)*comp_inter[[17]]+fac_inter[[17]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facd10.rda", interact= (-1)*comp_inter[[17]]+fac_inter[[17]])

g_metric_FacCompDenseSp10<-metrics_gjam(Rho_sign,comp=comp_inter[[17]], fac=fac_inter[[17]], only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDenseSp10$success_env,3)))
## Success rate for non-interacting: 0.343
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDenseSp10$success_comp,3)))
## Success rate for competition: 0.6
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDenseSp10$success_fac,3)))
## Success rate for facilitation: 0.667

HMSC

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd10.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[17]] + fac_inter[[17]])

h_metric_FacCompDenseSp10<-metrics_hmsc(Rho_sign,comp=comp_inter[[17]],fac=fac_inter[[17]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDenseSp10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDenseSp10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDenseSp10$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Competition +Facilitation 20 dense species

JSDM

HMSC

data<-sim_data$FacCompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd20.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[18]] + fac_inter[[18]])
h_metric_FacCompDenseSp20<-metrics_hmsc(Rho_sign,comp=comp_inter[[18]],fac=fac_inter[[18]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDenseSp20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDenseSp20$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDenseSp20$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Competition +Facilitation 5 sparse species

JSDM

Gjam

data<-sim_data$FacCompSparseSp5
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs5.rda",interact= (-1)*comp_inter[[19]]+ fac_inter[[19]] )
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facs5.rda", interact=  (-1)*comp_inter[[19]]+ fac_inter[[19]])

g_metric_FacCompSparseSp5<-metrics_gjam(Rho_sign,comp=comp_inter[[15]], only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparseSp5$success_env,3)))
## Success rate for non-interacting: 0.778
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparseSp5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparseSp5$success_comp,3)))
## Success rate for competition: 1

HMSC

data<-sim_data$FacCompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs5.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[19]]+ fac_inter[[19]])

FacCompSparseSp5<-metrics_hmsc(Rho_sign,comp=comp_inter[[19]],fac= fac_inter[[19]] ,only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(FacCompSparseSp5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(FacCompSparseSp5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(FacCompSparseSp5$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Competition +Facilitation 10 sparse species

JSDM

Gjam

data<-sim_data$FacCompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs10.rda",interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facs10.rda", interact=  (-1)*comp_inter[[20]] +fac_inter[[20]])

g_metric_FacCompSparseSp10<-metrics_gjam(Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparseSp10$success_env,3)))
## Success rate for non-interacting: 0.415
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparseSp10$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparseSp10$success_fac,3)))
## Success rate for facilitation: 0

HMSC

data<-sim_data$FacCompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs10.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[20]] +fac_inter[[20]])

h_metric_FacCompSparseSp10<-metrics_hmsc(Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparseSp10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparseSp10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparseSp10$success_fac,3)))
## Success rate for facilitation: 0

Environmental filtering + Competition +Facilitation 20 sparse species

JSDM

Gjam

data<-sim_data$FacCompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs20.rda",interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facs20.rda", interact= (-1)*comp_inter[[21]]+fac_inter[[21]])

g_metric_cmp_facs20<-metrics_gjam(Rho_sign,comp=comp_inter[[21]], fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_cmp_facs20$success_env,3)))
## Success rate for non-interacting: 0.363
cat(sprintf("Success rate for competition: %s\n", round(g_metric_cmp_facs20$success_comp,3)))
## Success rate for competition: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_cmp_facs20$success_fac,3)))
## Success rate for facilitation: 0.75

HMSC

data<-sim_data$FacCompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs20.rda" )
hm_conv(hm_mod)

Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[21]] +fac_inter[[21]])
h_metric_cmp_facs20<-metrics_hmsc(Rho_sign,comp=comp_inter[[21]],fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_cmp_facs20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_cmp_facs20$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_cmp_facs20$success_fac,3)))
## Success rate for facilitation: 0

Graph from paper

Visualization for true interactions

Facilitation

Competition Dense

Competition Dense

Competition + Facilitation Dense

Competition + Facilitation Sparse